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Abstract 

We present a new type of cluster algorithm that strongly reduces critical 
slowing down in simulations of vertex models. Since the clusters are 
closed paths of bonds, we call it the loop algorithm. The basic steps in 
constructing a cluster are the break-up and the freezing of vertices. We 
concentrate on the case of the F model, which is a subset of the 6-vertex 
model exhibiting a Kosterlitz-Thouless transition. The loop algorithm 
is also applicable to simulations of other vertex models and of one and 
two-dimensional quantum spin systems. 



PACS numbers: 02.70+d, 05.50+q, 75.10Jm, 68.35Rh 



Introduction 



Cluster algorithms [p], |2| are one of the few known ways to overcome critical slowing 
down in Monte Carlo simulations. Starting with l| and continuing with new ideas 
like H and ||, most of the successful algorithms have dealt with spin systems with 
two-spin interactions (see however ||). 

In vertex models |7], § the dynamical variables are localized on bonds, and the 
interaction is between all bonds meeting at a vertex. Furthermore there are con- 
straints on the possible bond variable values around a vertex. 

In this paper we present the loop algorithm, a new type of cluster algorithm 
applicable to vertex models. For usual spin systems most cluster algorithms start 
by "freezing" (also called "activating") or "deleting" bonds. Clusters are then sets 
of sites connected by frozen bonds. In the case of vertex models our idea is to define 
clusters as closed paths of bonds ("loops"). To construct such clusters, we have to 
perform operations at vertices that generalize the freeze-delete procedure. In this 
context we introduce the concept of break-up of a vertex. 

For the sake of clarity we concentrate on the F model, which is one of the simplest 
vertex models. We define it on an L x L square lattice. Vertices are located at lattice 
sites. The bond variables take the values ±1. They can be represented by arrows 
(e.g. +1 means arrow up or right, —1 means arrow down or left). At each vertex 
we have the constraint that the number of incoming arrows equals the number of 
outgoing arrows. Thus there are six different vertex configurations (six "vertices"), 
as shown in fig. 1. Their statistical weights w(i), i = 1, . . . , 6 are: 



c 



K i = l,2,3,4 



The coupling K > plays the role of inverse temperature. At K c = \a2 there is 
a Kosterlitz-Thouless transition. The correlation length is finite for K > K c and 
infinite for K<K C . 

In what follows we start by presenting our new loop algorithm. It turns out that 
there is one free parameter in the algorithm. We discuss how to choose an optimal 
value. Then we analyze the exponential autocorrelation times at K = K c and at 
K = K c /2. For the optimum algorithm we find a dynamical critical exponent of 
z(K c ) = 0.71(5) and z(K c /2) = 0.19(2). No critical slowing down is visible for the 
total energy. We briefly show how to generalize our algorithm to more general six 
and eight vertex models and how to use it for simulations of quantum spin systems. 



The Loop Algorithm 

If we regard the arrows on bonds as a vector field, the constraint at the vertices is 
a zero-divergence condition. Therefore every configuration change can be obtained 
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1 2 3 4 5 6 

Figure 1: The six vertex configurations. The labels 1, . . . , 6 follow standard 
conventions 



ul-lr 11-ur 

Figure 2: The two break-ups of a vertex: ul-lr (upper-left-lower-right) and 
11-ur (lower-left-upper-right) 



as a sequence of loop-flips. By "loop" we denote an oriented, closed, non-branching 
(but possibly self-intersecting) path of bonds, such that all arrows along the path 
point in the direction of the path. A loop-flip reverses the direction of all arrows 
along the loop. 

Our cluster algorithm performs precisely such operations, with appropriate prob- 
abilities. It constructs closed paths consisting of one or several loops without com- 
mon bonds. All loops in this path are flipped together. 

We shall construct the path iteratively, following the direction of the arrows. Let 
bond b be the latest addition to the path. The arrow on b points to a new vertex 
v. There are two outgoing arrows at v, and what we need is a unique prescription 
for continuing the path through v. This is provided by a break-up of the vertex v. 
In addition to the break-up, we have to allow for freezing of v. By choosing suitable 
probabilities for break-up and freezing we shall satisfy detailed balance. 

The break-up operation is defined by splitting v into two corners, as shown in 
fig. 2. At any corner one of the arrows points towards v , while the other one points 
away from v. Thus we will not allow e.g. the ul-lr break-up for a vertex in the 
configuration 3. A "corner flip" is a flip of both arrows. For a given break-up, 
we only allow the configuration changes resulting from independent corner flips. 
This preserves the zero divergence condition at v. Notice that a single corner flip 
transforms a vertex of weight 1 into a vertex of weight e~ K and vice-versa. Detailed 
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balance is satisfied with the following probabilities for choosing a given break-up: 
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P ul-lrW= * = 3 ' 4 > P \\-nA t ) = \ re ^ = 3,4 ; (2) 
[ r ? 

r is a free parameter for now. 

Freezing of a vertex means that its weight must not change. Since there are 
only two different vertex weights, we introduce two freezing probabilities. They are 
already determined by the requirement that for a given vertex configuration the sum 
of freezing and break-up probabilities must be one: 

P m _ J l-re* »= 1,2,3,4 . , 

^freeze W ~ j i_2r z = 5,6 ' 1 j 

The range of possible values for r is now obtained by requiring that all probabilities 
are between zero and one: 

< r < min(~,e~ K ) . (4) 

Assume now that we have broken or frozen all vertices. Starting from a bond bo, 
we proceed to construct a closed path by moving in the arrow direction. As we move 
from vertex to vertex, we always have a unique way to continue the path. If a vertex 
is broken, we enter and leave it along the same corner. If the vertex is frozen and 
of type 1, 2, 3 or 4, we pass through it on a straight line. At such vertices the path 
may be self-intersecting. Finally, if the latest bond b added to the cluster points to 
a frozen vertex v of type 5 or 6, the path continues both to the right and to the 
left of b, i.e. we start a new loop at v. The two loops have to be flipped together. 
In general, the zero-divergence condition guarantees that all loops will eventually 
close. 

The break-or-freeze decision for all vertices determines a unique partitioning of 
the lattice into closed paths that can be flipped independently. We choose to perform 
single cluster updates, i.e. we "grow" a single path from a random starting bond b , 
and flip it. The break-or-freeze decision is only needed for the vertices along the 
path. Thus the computer time for one path is proportional to the length of that 
path. 

It is easy to see that our algorithm is correct. The proof of detailed balance is 
completely analogous to that for other cluster algorithms (l|, || . The main ingredient 
here is that P n \_\ T and -Pvj_ ur already satisfy detailed balance locally. Furthermore, 
it is not difficult to see that any two allowed configurations can be connected by a 
finite number of cluster flips. Thus a finite power of the Markov matrix is ergodic. 

How do we choose an optimal value for the parameter r ? We have seen that 
freezing of a vertex of type 5 or 6 forces us to flip two loops together. If we had 
broken it up instead, we might have been allowed to flip the two loops independently. 
Thus more freezing leads to larger clusters. We conjecture that the least possible 
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freezing is optimal. This is confirmed by numerical tests (see below). From eq. (|3|) 
we then obtain 

r ° pt " \ e~ K K>K C • [b) 

By maximizing r we also minimize the freezing probability for vertices of type 1, 2, 
3 and 4. Notice that if we choose r = r opt , then for K <K C vertices of type 5 and 
6 are never frozen, so every path consists of a single loop. For K>K C on the other 
hand, vertices of type 1, 2, 3 and 4 are never frozen, so we do not continue a path 
along a straight line through any vertex. 

There are some distinct differences between our loop-clusters and more conven- 
tional spin-clusters. For spin-clusters, the elementary objects that can be flipped are 
spins; freezing binds them together into clusters. Our closed loops on the other hand 
may be viewed as a part of the boundary of spin-clusters (notice that the boundary 
of spin clusters may contain loops inside loops). It is reasonable to expect that in 
typical cases, building a loop-cluster will cost less work than for a spin-cluster. This 
is an intrinsic advantage of the loop algorithm. 

This can be exemplified nicely for the F model, where a spin-cluster algorithm 
- the VMR algorithm || - is also available. At K c one can see that if we use 
r = r o P t, loop-clusters are indeed parts of the boundary of VMR spin-clusters. Since 
flipping a loop-cluster is not the same as flipping a VMR cluster, we expect the two 
algorithms to have different performance. We found (see |J and the next section) 
that in units of clusters, the VMR algorithm is more efficient, but in work units, 
which are basically units of CPU time, the loop algorithm wins. At K c /2, where 
the loop-clusters are not related |T(| to the boundary of VMR clusters, we found 



the loop algorithm to be more efficient both in units of clusters and in work units, 
with a larger advantage in the latter. 



Performance 

We tested our new algorithm on L x L square lattices with periodic boundary condi- 
tions, both at the transition point K c and at \K C deep inside the massless phase. We 
carefully analyzed autocorrelation functions and determined the exponential auto- 
correlation time t. At infinite correlation length, critical slowing down is quantified 
by the relation [|TJ 

tocL z . (6) 

Local algorithms are slow, with z ~ 2. For comparison, we performed runs with 
a local algorithm that flips arrows around elementary plaquettes with Metropolis 
probability, and indeed found z = 2.2(2) at K — K c . 

In order to make sure that we do observe the slowest mode of the Markov matrix 
we measured a range of quantities and checked that they exhibit the same r. As in 
0, the slowest mode is strongly coupled to the sublattice energy. The two sublattice 
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L 


K = K C 


K = \K C 


8 


1.8(1) 


4.9(4) 


16 


3.0(2) 


5.6(2) 


32 


4.9(4) 


6.2(3) 


64 


7.2(7) 


7.4(3) 


128 


15.5(1.5) 


8.3(2) 


256 


20.5(2.0) 




z 


0.71(5) 


0.19(2) 



Table 1: Exponential autocorrelation time r at r = r opt , and the resulting 
dynamical critical exponent z. 

energies |J add up to the total energy. The constraints of the model cause them to 
be strongly anticorrelated. Within our precision the true value of r is not visible 
from autocorrelations of the total energy, which decay very quickly. Only for the 
largest lattices do we see a small hint of a long tail in the autocorrelations. A 
similar situation occurred in M, where, when decreasing the statistical errors, the 
decay governed by the true r eventually became visible. Note that as a consequence 
of this situation, the so-called "integrated autocorrelation time" KJ is much smaller 
than r, and it would be completely misleading to evaluate the algorithm based only 
on its values. 

We shall quote autocorrelation times r in units of "sweeps" [I]. We define a sweep 
such that on average each bond is updated once during a sweep. Thus, if r cl is the 
autocorrelation time in units of clusters, then r = r cl x < cluster size > /(2L 2 ). 
Each of our runs consisted of between 50000 and 200000 sweeps. Let us also define 

cl 

z c\ r ci ^ j^z ^ anc j a c i us t er s i ze exponent c by < cluster size > cx L c . We then 
have: 

z = z cl - (2 - c) . (7) 

Table 1 shows the autocorrelation time r for the optimal choice r = r opt . At 
K = \K C) deep inside the massless phase, critical slowing down is almost completely 
absent. A fit according to eq. || gives z = 0.19(2). The data are also consistent with 
z = and only logarithmic growth. For the cluster size exponent c we obtained 
c = 1.446(2), which points to the clusters being quite fractal. At the phase transition 
K = K c we obtained z = 0.71(5), which is still small. The clusters seem to be less 
fractal: c = 1.060(2). 

We noted above that at K = K c and for the optimal choice of r, the loop-clusters 
are related to the VMR spin-clusters. In |J we obtained for the VMR algorithm 
at K = K c the result z cl = 1.22(2), but we had c = 1.985(4), which left us with 
z — 1.20(2). In this case it is the smaller dimensionality of the clusters that make 
the loop algorithm more efficient. 
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K 


1 


/C 




0.500 


0.19(2) 




0.450 


1.90(5) 


1 TS 


0.400 


> 2.6(4) 


K c 


0.500 


0.71(5) 


K c 


0.475 


0.77(6) 


K c 


0.450 


0.99(6) 


K c 


0.400 


> 2.2(1) 



Table 2: Dependence of the dynamical critical exponent z on the parameter 
r. We use ">" where for our lattice sizes r increases faster than a 
power of L. 



As mentioned, no critical slowing down is visible for the integrated autocorre- 
lation time of the total energy. At K — K c , r int (E) is only 0.80(2) on the largest 
lattice, and we find z- mt (E) « 0.20(2). At K = \K a , r int (E) is 1.1(1) on all lattice 
sizes, so Zi n t(E) is zero. 

What happens for non-optimal values of r ? Table |2] shows our results on the 
dependence of z on r. z rapidly increases as r moves away from r opt . This effect 
seems to be stronger at \K C than at K c . We thus see that the optimal value of r 
indeed produces the best results, as conjectured from our principle of least possible 
freezing. 

In the massive phase close to K c , we expect ]10| that z(K c ) will determine the 



behaviour of r in a similar way as in ref. 0. To confirm this, a finite size scaling 
analysis of r is required. 



Generalizations and Outlook 

For the sake of clarity we have described our approach in terms of the F model 
only. It has however a much wider range of applicability. We will give a detailed 
description elsewhere ||10|| . Here we shall mention only a few highlights. 



Our "break-up of vertices" and subsequent path flip automatically satisfies the 
constraints of the F model. General six and eight vertex models 0j with arrow flip 
symmetry have related constraints. By using the framework of Kandel and Domany 
0] and the principle of minimal freezing, we can generalize the break-up operation 



TTj] to obtain efficient algorithms for these cases too. Algorithms for more general 
vertex models can be engineered along the same lines. 

Particularly promising is the possibility of accelerating Quantum Monte Carlo 
simulations ]lT|, [HJ. Quantum spin systems in one and two dimensions can be 
mapped into vertex models in 1 + 1 and 2+1 dimensions via the Trotter formula 
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and suitable splittings of the Hamiltonian [IT]. The simplest example is the spin 
\ xxz quantum chain, which is mapped into the 6-vertex model. For higher spins, 
more complicated vertex models result (e.g. 19- vertex model for spin one). 

For (2 + 1) dimensions, different splittings of the Hamiltonian can lead to geo- 
metrically quite different situations JT0|, [IT). We can e.g. choose between 6-vertex 
models on a complicated 2 + 1 dimensional lattice, and models on a bcc lattice with 8 
bonds (and a large number of configurations) per vertex. In all these cases, the con- 
straints are of a similar nature as in the F model, and our approach of constructing 
and updating clusters which are paths can be applied in a straightforward way. 

Notice also that in our approach it is easy to change global properties like the 
number of world lines or the winding number (see [IT|). 

Recently we received a paper by Wiese and Ying [|1^] on a different cluster 
algorithm for spin | quantum spin systems. After mapping to a vertex model similar 
to the one we refer to, they combine vertices into "block-spins" which are then 
used in a standard spin-cluster construction. This approach restricts the possible 
updates of the arrows. In our language, their clusters are sets of loops that are 
frozen together, i.e. that have to be flipped together. For some interesting cases, 
e.g. the one-dimensional Heisenberg ferromagnet and two-dimensional Heisenberg 
ferromagnet and anti- ferromagnet, the additional freezing leads to the problem of 
frustration for the block-spin clusters. We expect our algorithm to perform better 
in these cases, both because our clusters have less loops and because of the added 
flexibility offered by the possibility to optimize. 



Conclusions 

We have presented a new type of cluster algorithm. It flips closed paths of bonds 
in vertex models. Constraints are automatically satisfied. We have succesfully 
tested our algorithm for the F model and found remarkably small dynamical critical 
exponents. 

There are many promising and straightforward applications of our approach, 
to other vertex models, and to 1+1 and 2+1 dimensional quantum spin systems. 
Investigations of such systems are in progress. 
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